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A  technical  summary  of  this  project  is  contained  in  the  accompanying  presentation  materials. 
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of  propagation  of  atmospheric  light  strings  and  has  received  broad  exposure  amongst  academic, 
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propagation  experiments.  We  were  able  to  predict  the  phenomenon  of  femtosecond  light  string 
self-healing  on  collision  with  optically  opaque  scatterers.  A  fully  microscopic  theory  of 
incoherent  THz  emission  from  plasma  channels  left  in  the  wake  of  femtosecond  light  strings  has 
just  been  submitted  to  Physical  Review  Letters.  Our  UPPE  vector  ultra-short  pulse  propagator 
was  instrumental  in  identifying  a  new  class  of  nonlinear  X-waves  that  propagate  in  condensed 
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Sundstrom  ,  1972)  stability  theory  (also  known  as  Normal  Mode 
Analysis) 


AMR  object  coordinates  the  recursive  time-stepping  of  each  refinement  level 
Level  object  performs  memory  management  for  a  collection  of  subregions 
Interlevel  objects  handle  interpolation,  averaging  and  circulation  consistency 
Material  objects  encapsulate  an  algorithm  for  a  particular  material  model. 
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M.  Mansuripur,  A.R.  Zakharian  and  J.V.  Moloney,  "Transmission  of  light  through  small  elliptical  apertures  (Part  1)  ", 
Optics  and  Photonics  News,  March  Issue  pp44-49  (2004);  Part  2  April  issue  of  Optics  and  Photonics  News. 
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Self-healing  femtosecond  light  filaments 

M.  Kolesik1’2  and  J.V.  Moloney1’2 

1  Optical  Sciences  Center ,  University  of  Arizona  and 
2ACMS}  University  of  Arizona 

A  recent  experiment  (Courvoisier  et  al.,  Appl.  Phys.  Lett.  83(2003)213.)  indi¬ 
cated  that  filaments  created  in  femtosecond  high-power  pulses  propagating  in 
air  are  surprisingly  robust  when  interacting  with  microscopic  water  droplets. 
In  the  present  paper,  the  dynamics  of  the  filament-droplet  interaction  is  mod¬ 
eled  numerically.  Our  simulation  results  provide  further  insight  into  the  in¬ 
terplay  between  the  filament’s  core  and  the  wide  transverse  pedestal  of  the 
pulse.  It  is  shown  that  the  robustness  of  the  filament  comes  from  the  trans¬ 
verse  low-intensity  pedestal  that  controls  the  formation  of  the  central  hot  spot. 
Implications  for  penetration  of  wide,  high-power  beams  through  obscurants  are 
discussed. 

©  2003  Optical  Society  of  America 
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Filament  formation  in  ultrashort,  high-power  pulses 
continues  to  attract  interest  ever  since  its  first  experi¬ 
mental  observations  almost  a  decade  ago.1,2  While  there 
is  a  reasonably  good  understanding  of  the  processes  that 
control  the  single-filament  dynamics,3,4  the  very  high- 
power  regime,5  characterized  by  multiple  filaments  gen¬ 
erated  simultaneously,6  is  much  less  understood.  It  has 
been  predicted  by  computer  simulations  that  in  a  wide 
beam,  when  multiple  filaments  are  formed  simultane¬ 
ously,  the  low-intensity  background  of  the  pulse  plays 
a  crucial  role.7  It  serves  as  an  energy  reservoir  for  the 
emerging  and  decaying  filaments  and  thus  enables  the 
long-distance  propagation  of  the  resulting  “composite” 
pulse.  Clearly,  to  understand  the  dynamics  of  Tera-Watt 
scale  light  pulses,  that  can  easily  generated  hundreds 
of  filaments,  we  need  to  uncover  how  exactly  the  fila¬ 
ment  core  (center)  and  the  low-intensity  background  af¬ 
fect  each  other  during  propagation.  This  is  particularly 
important  for  propagation  through  obscurants,  such  as 
clouds. 

A  recent  experiment  by  Courvoisier  et  al.8  strongly  in¬ 
dicated  that  a  femtosecond  filament  is  extremely  robust 
when  it  interacts  with  microscopic  water  droplets.  They 
showed  that  filaments  with  sufficient  energy  can  survive 
head-on  collisions  with  droplets  as  large  as  95/mi  in  di¬ 
ameter.  This  ability  to  recover  from  the  collision  was 
attributed  to  energy  transfer  from  the  pulse’s  off-axis  re¬ 
gions.  Thus,  the  experiment  presents  another  strong  sup¬ 
port  for  the  spatial  replenishment  scenario  put  forward 
first  by  Mlejnek  et  al.3 

The  aim  of  this  letter  is  to  check  the  conclusions  of 
Ref.8  with  computer  simulations,  and  to  obtain  further 
insight  into  the  interaction  between  droplets  and  fila¬ 
ments.  Our  results  are  in  a  good  agreement  with  the 
experiment  by  Courvoisier  et  al.  Moreover,  we  show  that 
the  filament  robustness  has  little  to  do  with  its  alleged 
quasi-solitonic  character.  It  seems,  in  fact,  that  the  re¬ 
covery  after  the  collision  hardly  depends  on  the  filament 


structure  at  all.  Rather,  it  is  the  surrounding  low-energy 
background  of  the  pulse  that  replenishes  the  on-axis  re¬ 
gion  with  energy  and  continuously  builds  the  central  hot 
spot  quite  independently  of  what  happened  to  the  fila¬ 
ment  core  during  the  collision  with  the  droplet.  One 
can  say  that  the  short-scale  evolution  of  the  filament’s 
center  is  encoded,  or  “programed,”  in  the  low-intensity 
surrounding.  Our  simulations  also  indicate  that  the  sec¬ 
ondary  energy  losses  due  to  disturbed  waveform  in  the 
center  of  the  filament,  which  one  would  expect  to  result 
in  strong  radiation  of  the  energy  outward,  are  quite  small. 
It  seems  that  most  of  the  diffracting  field  is  re-captured 
within  the  pulse  and  continues  to  contribute  to  the  dy¬ 
namic  replenishment  in  the  later  propagation  stages. 

The  femtosecond  pulse  propagation  model  used  here  is 
based  on  a  simplified,  scalar  version  the  Unidirectional 
Pulse  Propagation  Equation.9  The  linear  air  susceptibil¬ 
ity,  needed  for  the  exact  linear  propagator,  is  constructed 
from  the  refractive  index  formula  by  Peck  and  Reed,10 
and  the  nonlinear  response  model  used  is  the  same  as 
in  Ref.3  Thus,  our  model  includes  effects  of  diffraction, 
self-focusing,  stimulated  Raman  scattering,  plasma  gen¬ 
eration  and  interaction  of  light  with  the  plasma.  As 
an  initial  condition,  we  use  a  loosely  collimated  (f=5m) 
Gaussian  pulse  (w  =  3mm)  of  140fs  duration  with  energy 
of  7m J.  The  unobstructed  pulse  creates  a  light  filament 
that  propagates  over  more  than  three  meters,  as  in  the 
experiment. 

The  experiment  indicates  that  only  the  geometric 
cross-section  of  the  droplet  plays  a  role,  while  its  op¬ 
tical  properties  are  unimportant.  Thus,  the  droplet  can 
be  modeled  simply  as  a  screen,  or  obstacle  of  appropriate 
size  in  the  simulation.  We  choose  our  model  droplet  to  be 
represented  by  a  screen  with  soft  edges.  Here  we  present 
results  for  the  screen  radius  of  50^m,  corresponding  to 
the  largest  droplets  studied  in.8 

One  intriguing  question  raised  by  the  experiment  is 
whether  or  not  it  matters  at  what  distance ,  after  the  fila- 
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Fig.  1.  Maximal  energy  fluence  as  a  function  of  the  propa¬ 
gation  distance  for  the  free  propagating  pulse  and  for  two 
pulses  hitting  droplets.  The  filament  dynamics  remains 
qualitatively  the  same  as  in  the  unperturbed  case.  For 
some  distance  after  the  droplets,  the  fluence  behaviors 
are  extremely  similar.  That  suggests  that  the  off-axis 
pulse  content  controls  the  central  hot  spot. 


mentation  onset,  the  filament  hits  the  droplet.  Some  re¬ 
searchers  in  the  field  believe  that  there  is  a  rather  delicate 
balance  between  self-focusing  and  de-focusing  forces  that 
holds  the  filament  together  and  produces  a  soliton-like 
propagating  structure.  We  believe  that  the  filament  goes 
through  different  stages  during  its  replenishment  cycles.3 
In  the  first  case,  the  filament  should  be  very  sensitive 
to  any  perturbations  in  its  central  region.  In  the  second 
case,  it  is  conceivable  that  droplets  encountered  at  differ¬ 
ent  stages  of  its  evolution  could  cause  different  degrees  of 
“damage”  to  the  filament  and  affect  differently  its  ability 
to  propagate.  We  have  therefore  studied  droplets  placed 
at  various  distances  chosen  to  intercept  different  dynamic 
states  of  the  pulse. 

A  quantity  that  is  relatively  easy  to  measure  in  exper¬ 
iments  is  the  fluence,  i.e.  the  total  energy  density  reach¬ 
ing  a  target  at  a  certain  spot.  We  compare  the  maximal 
fluences  recorded  over  several  meters  in  the  free  pulse 
and  in  pulses  that  encountered  droplets  at  z  =  0.4m  and 
z  =  1.0m  (Fig.  1).  From  the  point  of  view  of  the  filament 
dynamics,  these  are  two  qualitatively  different  locations. 
In  the  first  case,  the  filament  just  enters  its  first  replen¬ 
ishment  cycle,  and  has  a  much  simpler  waveform  than 
at  z  —  1.0m,  where  it  re- focuses  during  its  second  re¬ 
plenishment  cycle.  Figure  1  shows  that  in  both  cases  the 
filament  survives,  and  the  effect  of  the  droplet  on  the  fil¬ 
ament  propagation  distance  is  minimal.  Moreover,  just 
«  10cm  after  the  droplet,  there  is  a  striking  similarity 
between  the  fluence  behaviors  for  next  »  1m.  That  sug¬ 
gests  that  the  on-axis  pulse  shapes  depend  more  on  the 
content  of  the  off-axis  surroundings  than  on  the  central 
portion  of  the  filament  at  an  earlier  propagation  distance! 

Similar  to  maximal  energy  fluences,  the  plasma  genera¬ 
tion  is  only  slightly  affected  by  collision  with  the  droplet. 


Fig.  2.  The  overall  plasma  generation  remains  almost  un¬ 
affected  by  collisions  with  droplets,  independently  of  the 
pulse’s  evolution  stage  when  it  hits  the  droplet.  Under 
given  conditions,  the  plasma  generation  reflects  most  of 
the  energy  losses  of  the  pulse. 


This  is  demonstrated  in  Fig.  2:  The  total  number  of  free 
electrons  generated  per  unit  length  remains  practically 
the  same  in  perturbed  filaments.  Significant  differences 
only  occur  in  the  final  stages  of  the  propagation  where 
they  reflect  accumulated  deviations  between  their  gener¬ 
ating  optical  fields. 

Having  seen  the  robustness  of  the  filament  in  exper¬ 
imental  results,  and  corroborating  it  in  simulations,  we 
should  ask  the  question  of  how  much  the  droplet  actually 
“damages”  the  filament.  It  should  be  noted  that  our  large 
droplet  size  is  comparable  with  the  dimension  of  the  fila¬ 
ment’s  core,  and  the  collision  “annihilates”  a  significant 
portion  of  the  filament  itself.  However,  the  important 
point  is  that  the  total  energy  is  much  less  affected.  This 
is  shown  in  Fig.  3  where  we  plot  total  energies  of  the  free 
pulse  and  for  the  two  droplet-hitting  pulses.  Very  much 
like  in  the  experiment,  the  immediate  energy  loss  is  only 
several  percent  of  the  total. 

An  important  observation  here  is  that  the  subsequent 
energy  loss  rates  remain  almost  the  same  in  the  free  and 
perturbed  pulses.  If  there  was  a  significant  portion  of 
energy  radiated  out,  it  would  manifest  itself  as  an  energy 
decrease  compared  to  the  free  pulse.  (Note  that  we  use 
absorbing  boundary  to  keep  our  computational  domain 
reasonably  small).  Thus,  it  seems  that  the  diffraction 
rings,  or  ripples  created  by  the  obstacle  do  not  propa¬ 
gate  outside  of  the  beam,  but  are  captured  by  the  self- 
focusing  action  of  the  rest  of  the  pulse.  At  any  rate,  that 
portion  of  the  energy  that  is  radiated  out  is  probably  sig¬ 
nificantly  smaller  than  the  immediate  energy  loss  due  to 
the  collision  as  it  is  not  discernible  in  the  simulation. 

One  could  expect  that  as  a  result  of  taking  away  the  fil¬ 
ament’s  core,  the  self-focusing  force  that  holds  the  whole 
structure  together  diminishes,  and  diffraction  should  pre¬ 
vail.  However,  as  it  is  shown  in  Fig.  4,  the  hole  carved 
in  the  filament  core  by  the  droplet  closes  very  quickly. 
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Fig.  3.  Total  optical  field  energy  within  the  computa¬ 
tional  domain.  The  immediate  energy  lost  in  the  collision 
with  the  droplet  constitutes  only  several  percent  of  the 
total.  Overall,  energy  loss  rates  of  the  free-propagating 
and  perturbed  filaments  are  almost  the  same.  That  indi¬ 
cates  that  the  secondary  induced  energy  losses  (diffrac¬ 
tion)  are  negligible. 


Fig.  4.  Energy  fluence  with  and  without  the  droplet  at 
two  close  propagation  distances.  The  “hole”  in  the  pulse 
carved  out  by  the  droplet  is  shown  in  the  dotted  curve. 
Diffraction  together  with  the  absence  of  the  de-focusing 
effect  of  plasma  in  the  on-axis  region  causes  fast  “healing” 
of  the  hole.  Shortly  after  collision,  the  perturbed  filament 
exhibits  higher  on-axis  fluence  than  the  free-propagating 
filament.  Very  little  “ripple  effect”  is  visible. 

Just  10cm  after  the  droplet,  the  on-axis  concentration 
of  energy  is  stronger  than  in  the  the  free-propagating 
pulse.  This  fast  healing  of  the  hole  after  the  droplet  is 
due  to  the  combined  action  of  diffraction  and  expulsion 
of  the  light  from  the  off-axis  region  where  the  plasma  is 
still  generated.  Another  important  fact  is  that  there  is 
almost  no  “ripple  effect”  visible.  Comparing  transverse 
fluence  profiles  10cm  after  the  droplet,  it  is  difficult  to 
say  which  belongs  to  the  free-propagating  pulse.  This, 


we  speculate,  is  due  to  both  the  soft  edge  of  the  obstacle 
as  well  as  due  to  non-linear  capture  of  the  outgoing  light. 

The  above  observations,  though  obtained  only  for  a 
single  filament  in  a  radially  symmetric  simulation  have 
important  implications  for  high-power  pulses  with  mul¬ 
tiple  filaments  that  interact  with  obscurants.  The  fila¬ 
ments  can  overcome  collisions  with  droplets  due  solely  to 
the  low-intensity  background.  But  if  there  is  sufficient 
energy  in  the  latter,  then  such  characteristics  as  attained 
intensities,  energy  fluences,  and  plasma  generation  effi¬ 
ciency  are  largely  insensitive  to  perturbations.  Our  sim¬ 
ulations  even  indicate  that  the  secondary  energy  losses 
due  to  diffraction  are  small.  We  have  also  shown  that 
this  robustness  is  almost  entirely  provided  by  the  low- 
intensity  background,  while  the  structure  of  the  central 
core  of  the  filament  is  not  important.  These  findings  fur¬ 
ther  strengthens  the  dynamical  replenishment  scenario, 
and  together  with  the  fine  experiment  of  Courvoisier  et 
al.  exclude  any  purely  self-guiding  model. 

The  difficult  question  is  whether  the  filament  ro¬ 
bustness  implicates  robustness  of  the  whole  pulse  in 
the  multiple-filament  regime.  It  is  clear  that  the  low- 
intensity  background  cannot  be  treated  as  linear:  In  a 
wide  beam,  some  of  the  light  scattered  by  the  droplets 
may  be  recaptured.  It  will  be  important  to  under¬ 
stand  the  interplay  between  the  filaments  and  their  low- 
intensity  background. 

This  work  was  supported  by  Air  Force  Office  of  Scien¬ 
tific  Research  under  contract  F49620-03-1-0194.  Authors 
would  like  to  thank  Ewan  Wright  for  helpful  discussions. 
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Ultrashort  pulse  propagation  and  super  continuum  generation  in  tapered  and  microstructured 
optical  fibers  is  usually  simulated  using  the  corrected  nonlinear  Schrodinger  equation.  One  of  the 
underlying  approximations  is  that  of  wavelength-independent  effective  area  or,  equivalently,  constant 
nonlinear  coefficient  7.  In  very  thin  waveguide  structures  with  strong  light  confinement  the  validity 
of  such  an  approximation  comes  into  question.  In  this  paper  we  present  an  improved  model  in 
which  all  modal  properties  are  fully  taken  into  account  as  function  of  the  wavelength.  We  use 
comparative  numerical  simulation  to  identify  certain  regimes  in  which  an  improved  model  is  needed 
for  quantitatively  correct  results. 

PACS  numbers: 


I.  INTRODUCTION 

Recently,  novel  optical  fiber  structures  continue  to 
attract  much  interest.  In  particular,  photonic  crystal, 
micro-structured  [1,  2],  and  tapered  [3]  optical  fibers 
promise  a  wide  range  of  applications  and  are  therefore 
studied  intensively.  An  important  application  of  these 
waveguides  is  in  the  field  of  supercontinuum,  or  white- 
light  generation  [2-8]. 

Numerical  simulation  of  ultrashort  light  pulses  in 
micro-structured  and  tapered  fibers  is  usually  based  on 
the  corrected  Nonlinear  Schrodinger  Equation  (NLS) 
[6,  9,  10].  In  the  plethora  of  the  NLS  applications,  its 
corrected  form  used  in  fibers  is  probably  the  most  clean 
and  realistic  model.  Though  several  approximations  are 
necessary  to  derive  and  justify  this  equation,  in  the  con¬ 
text  of  optical  fibers  the  neglected  effects  are  usually  very 
small. 

The  corrected  NLS  has  been  also  used  to  gain  in¬ 
sights  into  the  dynamics  of  supercontinuum  generation 
in  micro-structured  and  tapered  fibers  with  diameters  as 
small  as  two  microns,  and  provided  results  in  good  qual¬ 
itative  agreement  with  experiments  [9,  10]. 

Recently,  we  pointed  out  that  even  thinner,  sub¬ 
micron  tapered  fibers  should  be  excellent  supercontin¬ 
uum  sources  [?  ].  However,  with  the  waveguide  diameter 
comparable  to  the  wavelength,  the  modal  electromag¬ 
netic  fields,  the  effective  mode  area  and,  consequently, 
the  nonlinear  coefficient  of  the  fiber  all  exhibit  an  ap¬ 
preciable  wavelength  dependence.  In  the  corrected  NLS, 
the  nonlinear  coefficient  is  taken  constant.  The  validity 
of  this  approximation  comes  into  question,  especially  for 
femtosecond  pulses  that  produce  supercontinuum  spec¬ 
tra  sometimes  spanning  more  than  an  octave  in  the  fre¬ 
quency  domain. 

It  is  the  goal  of  this  paper  to  clarify  whether  or  not  the 
corrected  NLS  can  be  reliably  applied  even  to  very  thin 
tapers  or  sub-micron  core  micro-structured  fibers.  We 
find  that  in  general  the  corrected  NLS  will  provide  results 
that  are  qualitatively  correct.  In  most  cases  that  the 
corrected  NLS  solutions  deviate  from  a  more  complete 


model  the  differences  are  too  small  to  be  significant  in 
the  qualitative  comparison  with  experimental  measure¬ 
ments.  However,  we  demonstrate  that  there  are  certain 
regimes  that  require  a  more  sophisticated  propagation 
model  for  semi-quantitative  simulations.  If  the  rapidly 
improving  experimental  techniques  for  measuring  the  su¬ 
percontinuum  pulses  allow  quantitative  comparison  with 
simulations,  and  we  believe  that  will  be  in  a  near  future, 
the  improved  propagation  model  will  be  required  in  many 
instances. 


II.  SUPERCONTINUUM  GENERATION  IN 
TAPERED  AND  MICRO- STRUCTURED  FIBERS 

Both  micro-structured  and  tapered  silica  fibers  have 
been  shown  to  be  promising  white-light  sources.  They 
allow  the  produce  super  continuum  spectra  often  span¬ 
ning  more  than  an  octave  at  relatively  low  powers.  How¬ 
ever,  optimizing  the  output  for  practical  applications  re¬ 
mains  a  problem  [7,  8].  In  general,  the  desired  spectral 
properties  are  large  width,  possibly  spanning  an  octave 
in  the  frequency  domain,  absence  of  deep,  wide  spectral 
gaps  and  smooth,  flat  spectral  intensity.  Simple  spectral 
phase  is  also  important  for  pulse  shaping.  Unfortunately, 
the  spectra  obtained  from  tapered  and  micro- structured 
fibers  do  not  usually  exhibit  all  these  properties.  Clearly, 
detailed  understanding  of  the  underlying  processes,  in¬ 
cluding  the  role  of  noise,  is  important. 

There  are  three  distinct  stages  of  the  supercontinuum 
generation  in  tapered  and  micro-structured  fibers.  These 
evolution  stages,  with  smooth  cross-overs  between  them, 
have  been  identified  in  several  computer  simulations,  in 
good  qualitative  agreement  with  experiments.  We  refer 
the  reader  to  Refs.  [6,  9,  10]  for  details,  and  summarize 
their  findings  together  with  our  simulational  observations 
next. 

First,  the  self-phase  modulation  is  the  prevailing  non¬ 
linear  effect  that  controls  the  initial  evolution  of  the 
pulse.  In  the  time  domain,  the  pulse  becomes  steeper, 
while  its  spectrum  broadens  in  a  nearly  symmetrical  fash- 
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ion.  During  this  initial  phase,  the  spectral  broadening  is 
relatively  slow  and  the  spectral  intensity  exhibits  charac¬ 
teristic  undulation  as  a  function  of  frequency. 

The  second  stage  is  characterized  by  fast  changes  in 
both  spectral  and  time  domains.  The  spectral  width 
grows  extremely  rapidly  over  propagation  distances  of 
the  order  of  millimeters  to  centimeters,  depending  on 
the  peak  pulse  power  and  the  effective  area  of  the  fiber. 
The  overall  spectral  shape  depends  on  concrete  condi¬ 
tions,  mainly  on  the  relative  position  of  the  central  pulse 
frequency  with  respect  to  the  chromatic  dispersion  pro¬ 
file  of  the  fiber.  However,  one  characteristic  trait  of  the 
spectrum  at  this  stage  is  that  it  has  relatively  simple 
structure,  exhibiting  several  broad  frequency  bands.  In 
certain  cases,  quite  flat  spectrum  without  deep  gaps  can 
be  achieved.  The  explosive  growth  of  the  spectral  width 
is  however  limited.  For  sufficiently  intense  pulses,  the 
total  attainable  with  is  rather  insensitive  to  pulse  peak 
power,  similarly  as  in  the  bulk  media;  The  spectral  with 
saturates  in  the  second  stage,  but  the  spectrum  continues 
to  evolve. 

During  the  third  stage,  a  complicated  fine-scale  spec¬ 
tral  structure  develops.  The  spectrum  gradually  evolves 
into  a  conglomerate  of  a  large  number  of  spectrally  nar¬ 
row  bands.  This  structure  continues  to  change,  with  its 
characteristic  scale  becoming  finer,  while  the  overall  “en¬ 
velope”  of  the  spectrum  remains  almost  unchanged.  Dur¬ 
ing  this  process,  the  coherence  is  gradually  deteriorating. 
It  is  this  stage  of  the  white-light  generation  process  that 
is  most  sensitive  with  respect  to  fluctuations  of  the  input 
pulse  and  other  noise  sources.  As  a  result,  the  detailed 
structure  of  the  supercontinuum  is  practically  random, 
completely  different  from  shot  to  shot.  After  sufficient 


propagation  length  in  the  fiber,  dispersion  effects  will 
eventually  take  over,  and  spectral  evolution  will  cease. 

Computer  simulation  have  been  providing  valuable  in¬ 
sights  into  the  dynamics  of  white-light  generation.  The 
computer  models  based  on  corrected  NLS  have  been 
found  in  a  very  good  qualitative  agreement  with  experi¬ 
mental  observations,  and  have  been  used  to  predict  some 
of  the  supercontinuum  properties.  However,  the  degree  of 
quantitative  accuracy  of  this  model  is  actually  not  quite 
certain  in  the  supercontinuum  generation  regime.  So  far 
the  great  sensitivity  with  respect  to  the  initial  condition, 
together  with  the  noise  amplification  and  limited  char¬ 
acterization  of  both,  the  input  and  the  output  prevent 
truly  quantitative  comparison  between  theory  and  exper¬ 
iment.  It  is  therefore  important  to  assess  the  accuracy  of 
the  model  also  by  other  means.  One  way  to  do  this  is  by 
comparison  with  a  more  complete  model,  and  that  is  the 
aim  of  the  next  section. 


III.  ULTRASHORT  PULSE  PROPAGATION 
MODELS 

With  the  detailed  derivation  provided  in  the  Appendix, 
here  we  only  describe  the  two  numerical  propagation 
models  employed  in  this  study.  Our  reference  model 
is  the  z-propagated  Unidirectional  Pulse  Propagation 
Equation  (UPPE,  Ref.  [11]  shows  derivation  of  the  time- 
propagated  version)  specialized  to  a  single  (transverse) 
mode  in  a  straight  silica  strand  of  radius  a  (see  Eq.  (25) 
in  Appendix): 
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Here,  the  electric  field  is  a  superposition  (in  the  frequency 
domain)  of  the  fundamental  modal  fields:  jE(r,t,  z)  — 
CUf(z)M(ricj)e~iu}t+i^^z .  This  sum  runs  over  a  dis¬ 
crete  set  of  angular  frequencies  u  that  correspond  to 
the  size  T  of  the  computational  domain  in  time.  ) 
is  the  fundamental  mode  propagation  “constant,”  and 
M(r,w)  is  the  normalized  vector  of  the  modal  elec¬ 
tric  field  calculated  exactly  for  all  cj’s.  Material  chro¬ 
matic  dispersion  properties  needed  for  that  calculation 
were  modeled  by  the  Sellmeyer  formula  [12]  for  the  sil¬ 
ica  index  of  refraction  with  neglected  losses.  We  use 
U2  =  2.7  x  10~16cm2/W  [13]  for  the  nonlinear  index, 
and  parameterize  the  response  function  of  the  stimulated 
Raman  effect  7 Z(r)  as  in  the  Intermediate  Broadening 


Model  [14].  For  notational  convenience  ,  we  include  in 
7 £(r)  also  the  instantaneous  optical  Kerr  effect  response. 

Equation  (1)  represents  a  large  system  of  ordinary  dif¬ 
ferential  equations  for  the  spectral  amplitudes  Cu,(z),  and 
various  standard  adaptive-step  ODE  can  be  used  to  solve 
it.  The  radial  integration  was  implemented  as  a  Gaus¬ 
sian  quadrature  over  radially  sampled  modal  fields.  Be¬ 
cause  of  the  frequency  dependence  of  the  modal  fields, 
this  integration  must  be  performed  for  each  right-hand- 
side  evaluation  and,  consequently,  this  simulation  is  an 
order  of  magnitude  slower  than  that  based  on  the  stan¬ 
dard  corrected  Nonlinear  Schodinger  Equation  which  we 
solve  in  the  spectral  domain: 


dzCu(z)  =  in 7  J  dte+^-W^Eit,  z)  jf°°  dTn{r)\E(t  -  r,  z) |2  ,  E(t, z)  =  J^  C»(z)e 
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We  term  this  equation  the  NLS  model  in  the  following 
comparison  with  the  reference  model.  It  is  essentially 
the  spectral  representation  of  the  corrected  Nonlinear 
Schrodinger  Equation.  The  only  difference  is  that  the 
NLS  is  usually  written  in  the  real-pace  representation 
with  the  free  propagator  (dispersion)  operator  most  of¬ 
ten  expressed  in  the  form  of  a  series  expansion  around 
the  reference  frequency.  The  reader  is  referred  to  Ap¬ 
pendix  for  details.  Here,  we  only  point  out  that  the  sim¬ 
plification  that  transforms  Eq.  (1)  into  Eq.  (2)  consists 
in  replacing  the  frequency-dependent  modal  fields  by  a 
single  modal  field  profile  at  a  chosen  reference  frequency 
wR:  M*{r,oJn).  It  is  best  to  choose  u;r  to 

be  the  carrier  frequency  of  the  incident  pulse,  of  course. 
Then,  with  the  radial  dependencies  of  E(r ,  t,  z)  fixed  this 
way,  one  can  perform  the  radial  integration  in  (1)  and 
pull  out  the  resulting  factor  outside  of  the  convolution 
integral  over  time. 

A  note  is  in  order  concerning  the  accuracy  of  our  ref¬ 
erence  model  itself.  Naturally,  although  in  certain  re¬ 
spects  better  than  the  NLS  model,  it  is  still  based  on  a 
number  of  simplifying  assumptions.  One  of  them  is  the 
neglected  frequency  dependence  of  the  nonlinear  index, 
which  should  have  a  similar,  although  smaller,  effect  as 
the  inclusion  of  the  frequency- dependent  modal  fields. 
Further,  the  finite  response  time  of  the  instantaneous 
Kerr  effect  response  may  affect  the  detailed  structure  of 
the  resulting  spectrum.  Also,  higher  order  modes  may  be 
excited  in  the  micro-structured  fibers  and  produce  extra 
broad  spectra.  Therefore,  while  our  comparison  will  help 
us  to  gauge  the  importance  of  a  particular  model  fea¬ 
ture,  it  may  not  tell  us  much  about  the  absolute  model 
accuracy. 

Finally,  it  is  to  be  noted  that  even  if  the  incident  field 
is  pure  single  (fundamental)  mode  (in  both  models),  the 
nonlinear  response  contains  terms  with  higher-order  an¬ 
gular  dependence,  and  in  principle  higher-order  modes 
are  always  generated.  However,  due  to  a  large  discrep¬ 
ancy  between  the  propagation  constants,  this  effect  is 
weak. 


IV.  NUMERICAL  SIMULATION  COMPARISON 

In  this  section  we  present  comparative  numerical  sim¬ 
ulation  simulations  for  the  two  models  described  above. 
Naturally,  it  is  impossible  to  make  generic  conclusions 
based  on  numerical  calculations  made  for  a  few  selected 
scenarios,  but  the  trends  we  observe  are  quite  clear,  and 
allow  us  to  make  some  useful  conclusions. 

Let  us  first  discuss  briefly  the  regimes  where  we  found  a 


very  good  agreement  between  the  reference  and  corrected 
NLS  models.  As  expected,  for  tapers  with  diameters  in 
the  vicinity  of  two  microns,  there  is  very  little  difference 
between  the  two.  The  overall  shape  of  the  supercontin¬ 
uum  spectra  are  so  similar  so  that  from  the  point  of  view 
of  comparison  with  current  experiments,  the  corrected 
NLS  model  is  fully  sufficient.  The  only  differences  in  the 
simulated  spectra  are  deviation  in  the  fine  spectral  struc¬ 
ture.  The  latter,  however,  is  to  a  large  extent  random, 
and  irreproducible  from  shot  to  shot.  It  is  also  sensitive 
to  noise  sources  included  in  simulation  and  as  such  only 
its  statistical  properties  could  be  in  principle  compared 
to  the  experimental  observations  of  multiple  single-shot 
spectra. 

On  the  other  hand,  for  sub-micron  diameter  silica 
strands,  the  two  models  deviate  from  each  other  signif¬ 
icantly  as  we  demonstrate  on  an  example  in  what  fol¬ 
lows.  This  is  not  to  say  that  the  standard  corrected  NLS 
should  be  automatically  replaced  by  the  more  complete 
reference  model  for  very  thin  fibers.  The  appropriate 
choice  of  model  will  depend  on  the  simulation  goals.  For 
qualitative,  exploratory  simulations,  the  corrected  NLS 
is  clearly  preferred  because  it  still  provides  qualitatively 
correct  picture  at  much  lower  computational  cost.  How¬ 
ever,  if  one  aims  at  more  accurate  results,  for  example 
when  studying  input-pulse  and/or  fiber  optimization  for 
supercontinuum  generation,  the  more  accurate  model  is 
recommended. 

Next,  we  show  results  for  simulation  of  the  100  fs  pulse 
centered  at  800  nm,  carrying  peak  power  of  10  kW,  prop¬ 
agating  in  a  short  length  of  a  straight  silica  strand  of 
radius  0.4/im.  Figure  1  shows  the  comparison  of  the  gen¬ 
erated  supercontinuum  spectra.  To  make  the  comparison 
easier  and  the  picture  more  readable,  we  have  smoothed 
both  spectra  (originally  normalized  to  maximum  spec¬ 
tral  intensity).  While  the  qualitative  picture  is  the  same 
in  both  models,  there  are  significant,  measurable  differ¬ 
ences,  too.  Most  importantly,  the  corrected  NLS  model 
exhibits  a  spectrum  narrower  by  tens  to  more  than  a 
hundred  nm  (depending  on  the  chosen  reference  level, 
of  course).  Also  the  spectral  gap  between  the  two  main 
spectral  peaks  appear  to  be  pronounced  more  in  the  NLS 
model.  Even  more  striking  differences  are  observed  in  the 
fine-scale  spectral  structure  shown  in  the  bottom  panel. 

Naturally,  the  difference  between  the  two  model  be¬ 
comes  even  more  pronounced  in  the  real  space.  Unlike  in 
the  spectral  space,  where  we  still  have  a  qualitative  simi¬ 
larity  (however,  let  us  keep  in  mind  that  the  comparison 
in  the  logarithmic  scale  tends  to  hide  differences),  the 
time-dependent  intensity  profiles  of  the  simulated  pulse 
deviate  quite  quickly.  This  is  shown  in  Fig.  2  for  two 
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wavelength  [nm] 


wavelength  [nm] 

FIG.  1:  Supercontinuum  spectra  generated  in  a  submicron- 
diameter  tapered  silica  strand.  Heavy  line:  corrected  NLS, 
thin  line:  reference  model.  Upper  panels  shows  smoothed 
spectra  in  order  to  reveal  the  tendency  of  the  deviation  be¬ 
tween  the  two  models  more  clearly.  The  lower  panel  shows 
the  detail  of  the  non-smoothed  spectra  revealing  completely 
different  fine-scale  structure. 


propagation  distances.  While  at  z  =  1.0cm,  one  can 
make  a  one-to-one  identification  of  the  pulse  strictures, 
at  z  =  2.0cm  the  time-profiles  appear  to  be  rather  differ¬ 
ent.  Of  course,  because  of  the  difference  in  the  spectral 
content,  the  deviation  in  the  real  space  will  continue  to 
grow. 


V.  CONCLUSION 

We  have  compared,  from  the  point  of  view  of  numerical 
simulation  of  the  ultrashort  pulse  propagation  in  tapered 
silica  waveguides,  the  widely  used  corrected  Nonlinear 
Schrodinger  equation  (NLS  model)  and  a  single-mode 
version  of  the  Unidirectional  Pulse  Propagation  Equation 
(reference  model)  with  emphasis  on  sub-micro  diameter 
silica  strands. 

While  the  NLS  model  always  provides  qualitatively 
correct  answers,  the  quantitative  deviations  from  the 
more  accurate  reference  model  become  significant  in  very 


local  pulse  time  [s] 


local  pulse  time  [s] 


FIG.  2:  Intensity  as  a  function  of  the  local  pulse  time  obtained 
from  the  standard  and  full  models.  Upper  panel  shows  the 
simulated  pulses  at  the  propagation  distance  z  —  1.0cm,  while 
the  lower  panel  is  for  z  =  2.0cm. 


thin  silica  strands.  Although  the  limited  number  of  sim¬ 
ulated  scenarios  prevent  us  from  drawing  general  conclu¬ 
sions  about  the  nature  of  the  deviations  between  the  two 
models,  we  believe  that  our  simulations  provide  useful 
indications  about  the  appropriate  usage  of  the  respective 
models.  Namely,  when  comparisons  with  an  experiment 
are  limited  to  multi-shot  spectra,  without  detailed  knowl¬ 
edge  about  the  input  pulse  as  it  is  usually  the  case  in  cur¬ 
rent  experiments  and  simulations,  the  currently  widely 
used  NLS  model  remains  the  method  of  choice.  It  pro¬ 
vided  qualitatively  correct  results  (as  compared  to  those 
obtained  from  the  reference  model)  in  all  simulations  we 
have  performed,  including  the  extreme  cases  of  very  thin 
tapers  and  high  peak  powers.  The  low  computational 
cost  of  the  NLS  model  makes  it  useful  for  exploratory 
runs  even  in  the  regimes  that  require  the  improved  model. 

Our  simulations  show  that  the  two  models  provide 
significantly  different  results  for  ultrashort  pulse  prop¬ 
agation  in  very  thin,  sub-micron  diameter  tapered  fibers 
when  the  spectral  broadening  causes  significant  portion 
of  the  pulse  energy  to  be  transferred  toward  the  long- 
wavelength  end  of  the  spectrum.  Namely  at  long  wave¬ 
lengths,  the  frequency  dependence  of  the  modal  fields 
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is  the  strongest  making  the  notion  of  the  frequency  in¬ 
dependent  nonlinear  coefficient  7  less  accurate.  Thus,  in 
cases  when  accurate  results  are  sought,  for  example  when 
using  simulations  to  optimize  the  input  pulse  and/or  the 
properties  of  the  fiber  for  supercontinuum  generation,  the 
more  accurate  model  should  be  used.  Moreover,  we  be¬ 
lieve  that  in  the  near  future  the  experimental  outputs 
in  the  field  of  SC  generation  in  tapered  and  microstruc- 
tured  fibers  will  become  suitable  for  detailed  comparison 
with  theory,  and  thus  will  make  it  possible  to  quantify 
the  accuracy  of  various  theoretical  models. 


VI.  APPENDIX 


The  sum  runs  over  all  transverse  modes  and  a  discrete  set 
of  angular  frequencies,  the  latter  corresponding  to  a  fi¬ 
nite,  large  normalization  “volume”  T  in  the  time  domain. 
In  what  follows,  we  use  the  shorthand 


for  all  time-domain  integrations  unless  integration 
bounds  are  shown  explicitly.  To  keep  notation  short,  we 
also  use  the  convention  that  if  modal  fields  are  used  with¬ 
out  showing  their  explicit  arguments,  the  time-dependent 
and  propagation  phase  factors  are  understood  to  be  ab¬ 
sorbed  into  modal  fields: 


For  the  sake  of  completeness,  as  well  as  for  reader’s 
convenience,  we  provide  below  a  detailed  derivation  of  the 
propagation  equations  for  femtosecond  pulses  in  straight 
silica  stands.  First  we  follow  a  textbook  derivation  of  a 
general  propagation  equation  for  the  pulse  optical  field 
expressed  as  a  superposition  of  guided  modes. 

Consider  a  nonmagnetic,  dispersive  medium  with  the 
relative  permitivity  e  that  only  depends  on  the  transverse 
coordinates  x,y  and  on  the  angular  frequency  u 

e  =  e((j,x,y)  ,  fi  =  no  ■  (3) 

Nonlinearity  and  all  other  effects  that  we  can’t  include 
exactly  will  be  lumped  in  the  polarization  P  in  the  ma¬ 
terial  costitutive  relation: 


Sm,u>  =  €mAx>y)eiMW)Z~iU,t 

i,  =  (7) 

The  orthogonality  relation 

J  z-[£mtw  "Hm,w  dxdy  =  26m}nNmilj  (8) 

will  be  used  below. 

Having  fixed  the  notation  we  can  derivation  the  z- 
propagated  Unidirectional  Pulse  Propagation  Equation. 
What  follows  is  a  textbook  method  based  on  the  reci¬ 
procity  relation.  We  start  from  Maxwell’s  equations 

dtP  -\-  £o  *  E  =  V  x  H 

-HodtH  =  Vxl  (9) 


D  =  €06  *  E  +  P  .  (4) 

Let  our  initial  data  be  given  in  the  x,y,t  domain. 
Then,  the  electromagnetic  field  of  a  pulse  propagating 
along  the  2~axis  can  be  expressed  as  a  superposition  of 
the  (fiber)  modal  fields 


which  we  scalar-multiply  by  the  complex  conjugate 
modal  fields,  including  their  time  and  2  phase  factors 

£m,u>'dtP  “keo SmiUrdt€  *  E  =  4.u,V  X  H 

=t,w.Vx£  (10) 

Let  us  re-arrange  both  right-hand  sides  as: 


e^.dtP  +  e0£^.dte  *  E  =  V. [if  x  4,J  +  H.[V  x  4,J 

-Ho =  V.[S  x  7C,J  +  E.[V  x  7C,J  (11) 

and  use  Maxwell’s  equations  for  the  conjugate  modes  that  appear  in  the  last  terms: 

4 ,„.dtP  +  eo4,a,A*  *  E  =  V.[S  x  4,J  - 

-HoH*m,u.dtH  =  V.[£  x  74, J  +  e0E.dte  *  4>w  (12) 

Next,  subtract  the  two  equations  and  integrate  over  the  whole  xyt  domain,  supposing  that  fields  vanish  at  infinity: 

J  S^.dtPdxdydt  =  dz  J  z.[H  x  £^J\dxdydt  -  dz  J  z.[E  x  H*m  u}dxdydt  (13) 

Finally,  keeping  in  mind  the  implicit  2,  t  dependence  of  the  modal  field,  insert  the  modal  expansion  (5)  and  use  the 
orthogonality  relation  (8)  to  obtain  the  evolution  equation  for  the  expansion  coefficients. 

dzAmtU  =  j  ,(*,„).£(*,  y,t)  dxdydt 


(14) 
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Next  we  specialize  the  propagation  equation  for  a  silica  strand  with  a  diameter  2 a.  Neglecting  the  nonlinear  interac¬ 
tions  outside  of  the  taper,  as  well  as  all  higher-order  modes,  the  propagation  equation  reads 

dzAw{z)  =  j  rdr  Jo  d</>  J  dtP*(r,<t>,t,z).P(r,4>,t,z)  (15) 

[ 


where  the  fundamental  mode  field  is 

f£r\  (£r{r,w)  COs(<j>)\ 

\  £</,  =  (  £<fi(r,u>)  sin(</>)  exp (~iu>t +  i/3(u)z)  (16) 

\£j  \£z(r,u>)cos(<f>)/ 

In  silica,  for  the  intensities  typically  encountered  in 
fibers,  the  polarization 

P(r,  <t>,  t,  z )  =  e0  A x(r,  <f>,  t,  z)  E{r,  <f>,  t,  z)  (17) 
is  expressed  through  the  local  change  of  the  susceptibility 

A*(r,  4>,  t,  z)  =  2 nhn2  [  drE2(r,  <f>,t-T,  z) TZ(r)  (18) 
Jo 

where  7£(r)  stands  for  the  normalized  Kerr/Raman  re¬ 
sponse  memory  function.  Since  we  neglect  coupling  to 
other  than  fundamental  modes,  we  need  to  consider  only 
the  axially  symmetric  component  of  Ax(r,  <£,£>*)  which 


is 

i  r°° 

A x(r,  t,  z)  =  2 nbn2-  /  dr[  V  E2(r,t-  r,  z)]K(t) 

2J0  c=r,*,z 

(19) 

We  have  used  the  fact  that  E^(r,  <j>,t  -  r ,  z)  is  a  super¬ 
position  of  fundamental  modes  (at  different  frequencies) 
and  thus  its  angular  dependence  is  proportional  to  ei¬ 
ther  cos(0)2  or  sin (</>)2,  depending  on  the  component  c. 
Therefore,  its  ^-independent  component  is  one  half. 

Thus,  in  the  single-mode  approximation,  the  polariza¬ 
tion  azimutal  dependence  is  the  same  as  of  the  one  of  the 
fundamental  mode: 

P(r,(j),t,z)  =  e0  Ax(r,t,z)  '^Aw{z)£w(r,  <f>,t,z)  (20) 

XV 

When  we  insert  this  expression  in  Eq.  (15),  the  angular 
integration  will  yield 


dzAw(z) 


iujireo 

2NU 


dte+iut-il3^z£*(r)E{r, 


t,  z)Ax(r,t,z) 


r 


(21) 


For  numerical  solver  implementation  it  is  convenient 
to  use  the  normalized  modal  fields.  Together  with  the 
corresponding  modal  propagation  constant,  the  modal 
fields  are  pre-calculated  exactly  for  a  range  of  angular 
frequencies  io\ 

P(aj)  ,  Mu(r)  =  £u{r)l^/K  (22) 

Similarly,  the  expansion  coefficients  are  replaced  by  the 
normalized  ones: 

Cu{z)  =  AM)VK,  (23) 


and  the  fields  expansion  then  reads 


2?(r,  t,z)  =  J2  Cu(z)Mu(r)e-iut+m“)z  (24) 


In  the  above  formulas,  vectors  are  expressed  in  the 
cylindrical  coordinates  and  the  angular  dependencies  are 
omitted.  With  this  normalization,  the  propagation  equa¬ 
tion  reads 


dzCu{z) 


%U)  7T6o 
2 


dte+^'-W^Ml  (r)  J(r,  t,  z)AX(r,  t,  z) 


(25) 


This  is  our  reference  model.  At  each  integration  step,  the  radial  integration  is  performed  using  Gaussian  or  Lo- 
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bato’s  numerical  scheme.  Thus,  the  evolving  pulse  field 
must  be  sampled  at  multiple  radial  locations.  Just  eight 
samples  turn  out  to  be  sufficient,  because  the  radial  field 
profiles  are  smooth.  Nevertheless,  the  radial  integration 
makes  this  method  order  of  magnitude  slower  than  the 
traditional  one-dimensional  corrected  NLS  which  we  de¬ 
rive  next. 

The  additional  approximation  we  need  to  make  to 
transform  the  above  propagation  equation  into  NLS-type 
equation  is  that  of  frequency  independent  modal  fields. 
Note,  that  the  modal  fields  enter  the  nonlinear  coupling 


term  at  four  generally  different  frequencies.  However, 
we  replace  (r)  — ►  A4^(r)  for  every  w,  with  fi  stand¬ 
ing  for  a  chosen  reference  angular  frequency.  Then,  the 
radial  dependency  can  be  isolated  from  the  field’s  time 
dependence.  Factoring  out  the  nonlinear  coefficient 

7  =  7rf°^t’n2-  rdr(M*n(r).Mci(r))1 2 3  (26) 
leads  to  the  corrected  NLS  written  in  the  spectral  domain 


dzCu(z)  =  iu 7  J  d te+^-W^’Eit, z)  jf°  ArTl{T)\E{t  -  r,  z) |2  ,  E(t,  z)  =  £  Cu (2)e-<“‘+^H  (27) 


This  we  refer  to  as  the  corrected  NLS  and  use  in  our 
comparative  numerical  simulations.  To  show  its  equiva¬ 
lence  with  the  usual  form  written  in  the  real-space  do¬ 
main,  one  needs  to  make  a  few  further  steps.  First, 
the  expansion  coefficient  are  redefined  to  absorb  the  free 
propagation  phase  factors  and  we  switch  to  the  envelope 
representation: 

E(t,  z)  =  A(t,z)e~iuot  = 

UJ 

(28) 

Second,  we  transform  the  evolution  equation  into  real- 

space  domain.  This  is  done  by  applying 

to  each  term  in  the  equation.  In  the  process,  (3(u> )  is 


interpreted  in  terms  of  Taylor  expansion  around  u;0,  and 
its  argument  is  replaced  by  u  u>o  4*  idt.  Then  the 
resulting  equation  reads 

dzA{t,  z)  =  i(3( vo  +  idt)A(t ,  z)  +  fy(u/o  +  idt)A7l  *  \A\2 

(29) 

Here,  /5(cj0 -H^t)  is  a  shorthand  for  a  differential  operator 
formed  from  the  Taylor  expansion  of  /5(u).  Finally,  to  get 
rid  of  the  first  two  terms  in  the  /5-expansion,  we  factor 
out  and  change  to  the  coordinate  frame  moving 

with  the  group  velocity  vg  =  l/(0w/?(o;))  evaluated  at 
u>0:  A(t  -  z/vg,z)  =  A{t,  which  gives  us  the 

corrected  NLS  equation  for  the  pulse  envelope 


dzA{t,z)  =  Yi\^P^(idt)nA{t,z)+iu>Ql{l+iuJ^dt)AK*\A\2  (30) 

1  n!  uu >n 

n=2 


r 


Note  that  u^o  is  an  “arbitrary”  parameter  introduced  in 
the  derivation,  and  recall  that  the  nonlinear  coefficient 
7  also  contains  similar  “degree  of  freedom,”  namely  f2. 
Naturally,  both  are  usually  chosen  equal  to  the  central 
frequency  of  the  initial  pulse. 

Thus,  equation  (27)  is  equivalent  to  the  standard  evo¬ 


lution  equation.  However,  the  latter  doesn’t  offer  any 
advantage  from  the  numerical  point  of  view.  Namely,  it 
requires  to  approximate  the  dispersion  relation  /5(w)  by 
a  polynomial,  and  the  higher  order  derivatives  are  best 
treated  in  the  spectral  domain  anyway.  It  is  therefore 
simpler  to  use  equation  (27). 
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Supercontinuum  generation  in  sub-micron  diameter  tapered  fibers 
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Super  continuum  (SC)  generation  in  tapered  fibers  with  taper  diameters  in  the 
sub-micron  range  is  investigated  numerically.  Such  thin  fibers  exhibit  two 
zero  dispersion  wavelengths  (ZDWs)  in  the  visible  range,  and  we  propose  and 
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improved  SC  generation. 
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Tapered1  and  microstructured2’3  optical  fibers  are 
promising  media  for  supercontinuum  (SC)  generation, 
as  wide  optical  spectra  can  be  generated  at  relatively 
low  powers  in  short  lengths  of  these  fibers. 1,4-6  How¬ 
ever,  the  quality  of  the  generated  white  light  remains 
a  problem.7,8  Ideally,  a  SC  source  should  provide  a  wide, 
flat  spectrum,  without  fine-scale  structures  or  spectral 
gaps.  For  generation  of  frequency  combs9,10  the  spec¬ 
trum  must  span  a  whole  octave  of  frequencies,  and  for 
pulse  compression  and  shaping,  the  coherence  properties 
and  the  phase  structure  of  the  supercontinua  are  impor¬ 
tant.6  Unfortunately,  tapered  fibers  and  microstructured 
fibers  normally  produce  continua  that  lack  most  of  these 
desired  properties.  Choosing  the  optimal  fiber  and  opti¬ 
mization  of  the  input  pulse  are  therefore  important.6,11 

One  of  the  attractive  features  of  the  tapered  silica  fiber 
is  that  it  exhibits  zero  dispersion  wavelength  (ZDW)  in 
the  visible  region  for  diameters  around  2  ^m.  Experi¬ 
ments  and  computer  simulations  have  therefore  concen¬ 
trated  on  tapers  with  dimensions  in  this  range.  The  aim 
of  this  Letter  is  to  demonstrate,  through  computer  sim¬ 
ulation,  that  tapered  fibers  with  sub-micron  diameters 
have  properties  that  are  well  suited  for  the  generation  of 
high-quality  SC  spectra,  in  particular,  they  exhibit  two 
ZDWs.  This  feature,  together  with  smaller  effective  ar¬ 
eas,  can  be  utilized  in  optimizing  the  tapered  fiber  for 
relatively  smooth  SC  spectra  .  Supercontinuum  gener¬ 


ation  around  the  second  ZDW  has  been  investigated  in 
the  infra-red  region.4  We  show  that  the  quality  of  the 
SC  generated  from  a  pulse  launched  in  the  vicinity  of  the 
second  ZDW  in  the  visible  frequency  region  is  superior 
to  that  obtained  from  a  thicker  fiber  utilizing  the  first 
ZDW.  Smooth  spectra  spanning  almost  an  octave  over 
the  visible  region  without  deep  gaps  can  be  obtained  in 
much  shorter  lengths  of  the  tapered  fiber. 

Another  potentially  useful  regime  is  realized  by  plac¬ 
ing  the  input  pulse  spectrum  symmetrically  between  two 
ZDWs.  The  symmetric  chromatic  dispersion  “landscape” 
then  gives  rise  to  a  spectrum  with  the  energy  concen¬ 
trated  symmetrically  at  the  extreme  low  and  high  fre¬ 
quencies.  Such  spectra  can  be  generated  at  low  energies 
and  thus  should  be  usable  for  frequency  combs10. 

In  the  present  qualitative  study,  we  concentrate  on  the 
tapered  section  of  the  fiber  only,  and  model  it  as  a  sim¬ 
ple  silica  strand  of  radius  a.  To  obtain  an  equation  that 
describes  the  evolution  of  the  pulse  along  the  strand‘s 
2-axis,  one  can  start  from  a  textbook  formula  for  mode 
excitation12  with  the  nonlinear  polarization  in  the  role 
of  the  “source.”  We  neglect  all  modes  except  the  funda¬ 
mental  one,  and  retain  only  the  symmetric  component  of 
the  nonlinear  response.  After  integrating  out  the  (funda¬ 
mental  mode)  polar-angle  dependence,  we  arrive  at  the 
following  propagation  equation  written  in  the  frequency 
domain 


dzCu(z)  =  /0  rdrM*{r,  w)e-tf  M*  t+°°  d tE(r,  t,  z)eiut  /°°  drn6n2ft(r)£2(r ,  t  -r,z)  .  (1) 
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Here  E(r,t,z)  =  f  dvCb,(z)/Z(r,v)e-ib*+Wb’)x  is  the 
electric  field  (with  the  polar-angle  dependence  sup¬ 
pressed),  is  the  frequency-dependent  propagation 

constant,  and  M(r,u)  is  the  normalized  vector  modal 
field  at  angular  frequency  oj.  These  quantities  are  cal¬ 
culated  exactly  for  a  silica  strand  over  a  sufficiently 


wide  range  of  frequencies.  We  use  a  Sellmeier  for¬ 
mula13  for  the  silica  index  of  refraction,  nonlinear  index 
ri2  =  2.7  x  10"16cm2/W.14  1Z(r)  stands  for  the  response 
function  of  the  stimulated  Raman  effect  (parameterized 
as  in  the  Intermediate  Broadening  Model  of  Ref.  [14]) 
that  includes  also  the  instantaneous  Kerr  effect  contri- 


2 


button.  These  coupled  equations  for  the  spectral  ampli¬ 
tudes  (normalized  to  unit  power)  are  solved  using  a 
standard  ODE  solver  routine. 

In  this  work  we  use  the  approximation  that  the  ex¬ 
act  frequency-dependent  modal  fields  can  be  replaced 
by  those  evaluated  at  the  center  frequency  of  the  input 
pulse  jM(r,  a;)  — >  With  this  simplification,  our 

model  reduces  essentially  to  the  widely  used  corrected 
Nonlinear  Schrodinger  Equation,15”17  the  only  difference 
being  that  the  nonlinear  coefficient  7  is  calculated  from 
the  vector  modal  field,  i.e.  not  from  the  “scalar”  effective 
area.  It  will  be  shown  in  a  separate  publication  that  for 
the  purpose  of  qualitative  studies  like  the  present  one, 
this  model  is  fully  sufficient.  Note  however,  that  in  the 
very  thin  tapers  there  are  certain  regimes  in  which  use  of 
the  full  equation  with  spectrally  resolved  modal  fields  is 
necessary  to  obtain  quantitatively  correct  results. 

When  the  SC  is  generated  in  a  tapered  or  a  mi- 
crostructured  fiber  there  are  three  distinct  stages  of 
propagation;6’16*17  In  the  first  stage,  self-phase  modu¬ 
lation  (SPM)  dominates  the  propagation  and  broadens 
the  spectrum  in  a  nearly  symmetrical  way.  In  the  second 
stage  at  larger  distances,  explosive  spectral  broadening 
occurs  and  the  spectral  width  of  the  SC  generation  satu¬ 
rates.  In  the  third  stage,  fine-scale  structures15,18  develop 
in  the  spectrum  and  coherence  is  gradually  lost  with  in¬ 
creasing  propagation  distance.7,8,16  To  obtain  wide  and 
coherent  spectra,  it  is  therefore  necessary  to  choose  an 
optimal  propagation  length  prior  to  the  third  propaga¬ 
tion  stage,6  and  we  have  chosen  the  propagation  lengths 
appropriately  here. 

To  proceed  we  compare  SC  generation  in  two  tapered 
fibers,  the  first  being  a  “standard,”  2.4/xm  diameter  ta¬ 
per  that  has  ZDW  located  close  to  A  =  780nm,  and  the 
second  a  thin  taper  with  a  radius  of  just  a  =  0.325^m. 
The  GVD  of  these  two  tapers  is  depicted  in  Fig.  1.  Note 
that  the  first  (and  only)  ZDW  of  the  thick  fiber  nearly 
coincides  with  the  second  ZDW  of  the  thin  fiber,  the  sec¬ 
ond  ZDW  vanishing  as  the  taper  diameter  is  increased. 
A  50fs  duration,  6kW  peak  power  pulse  is  launched  at 
A  =  780nm  in  both  fibers.  We  use  a  hyperbolic  secant, 
unchirped  pulse  for  this  comparative  simulations.  Natu¬ 
rally,  the  real  input  pulse  will  be  more  complex  than  that, 
but  for  a  more  realistic  comparison  of  SC  generation  in 
the  two  tapers,  one  should  include  the  input  transient 
taper  region  and  compare  results  for  inputs  optimized  for 
each  of  the  fibers.  Such  an  optimization  should  include 
shaping  of  the  pulse  and  is  much  more  complex  task  not 
addressed  in  this  work. 

Figure  2  compares  the  “optimal”  spectra  obtained  for 
the  two  fibers,  where  the  optimization  consists  of  choos¬ 
ing  the  interaction  length6  just  after  the  spectral  width 
saturates  and  before  it  develops  fine-scale  structure  that 
would  indicate  gradual  loss  of  coherence.  Further  propa¬ 
gation  in  the  taper  would  not  increase  the  useful  spectral 
width,  but  would  degrade  the  quality  of  the  spectra. 

The  key  observation  from  Fig.  2  is  that  the  spectrum  gen¬ 
erated  in  the  thicker  silica  taper  is  narrower  and  exhibits 


wavlength  [micron] 


Fig.  1.  Group  velocity  dispersion  of  the  fundamen¬ 
tal  mode  in  silica  strands  of  radius  a  —  0.325/im  and 
a=  1.2  pm. 


Fig.  2.  Supercontinuum  spectra  generated  in  two  dif¬ 
ferent  tapered  fibers.  Full  line:  0.325/mi  radius  taper, 
propagation  distance  z  —  1cm.  Dashed  line:  1.2/im  ra¬ 
dius  taper,  propagation  distance  z  =  20cm.  Dotted  line 
shows  the  input  spectrum. 


much  more  of  the  unwanted  spectral  structure  (dashed 
line),  in  comparison  to  the  spectrum  from  the  thinner 
fiber  (solid  line).  In  particular,  for  the  thicker  taper  the 
spectral  intensity  variations  are  deeper  and  the  gaps  in 
the  spectrum  are  wider  than  in  the  spectrum  obtained 
from  the  thin  strand.  The  greater  width  and  better  qual¬ 
ity  of  the  SC  obtained  from  the  thinner  strand  is  most 
likely  mainly  due  to  the  effectively  stronger  nonlinearity 
given  by  the  smaller  radius.  The  nonlinear  length  is  sig¬ 
nificantly  shorter  and  the  spectrum  is  generated  faster, 
which  gives  less  time  to  develop  the  fine-scale  structure 
and  spectral  gaps. 

The  computed  SC  spectrum  from  the  thin  taper  in 
Fig.  2  doesn’t  quite  span  a  whole  octave,  but  it  is  con¬ 
ceivable  that  the  result  can  be  improved  with  a  suitably 
shaped  input  pulse.  Our  preliminary  attempts  to  opti- 
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Fig.  3.  Supercontinuum  generation  from  a  pulse  launched 
in  the  middle  between  two  zero-dispersion  wavelengths 
in  a  a  —  0.25/mi  radius  taper  after  propagating  for  z  = 
lcm.  Further  propagation  leads  to  flattening  of  the  high- 
and  low-frequency  segments  of  the  spectrum,  but  fine- 
scale  spectral  undulations  and  structure  starts  to  develop 
quickly. 

mize  the  input  pulse  by  tuning  the  chirp  only  showed 
that  while  this  can  increase  the  spectral  width,  it  also 
results  in  a  spectrum  less  flat  with  deeper  gaps.  Thus, 
an  efficient  optimization  will  need  to  include  more  com¬ 
plex  pulse  shaping.  As  pointed  out  the  simulation  then 
must  include  the  whole  propagation  path  of  the  pulse 
with  changing  nonlinear  interactions  in  the  tapered  re¬ 
gions. 

The  result  we  want  to  emphasize  is  that  our  compara¬ 
tive  simulations  strongly  suggest  that  sub-micron  diam¬ 
eter  silica  strands  will  most  likely  be  better  sources  of 
SC  spectra  spanning  the  whole  visible  region  than  the 
d  «  2pm  diameter  fibers  used  in  number  of  recent  exper¬ 
iments. 

Another  possibility  for  visible  SC  generation  in  sub¬ 
micron  silica  strands  is  demonstrated  in  Fig.  3.  The  ini¬ 
tial  pulse  at  A  =  520nm  is  launched  in  the  middle  between 
the  two  ZDWs  of  a  a  =  0.30/im  radius  silica  strand.  The 
dispersion  profile  is  roughly  symmetric  around  this  cen¬ 
tral  wavelength,  and  the  value  of  the  dispersion  over  a 
wide  frequency  range  is  still  quite  small.  The  result  is  a 
highly  symmetric  spectrum  with  most  of  its  energy  con¬ 
centrated  at  low  and  high  frequencies.  The  generated 
SC  spectrum  spans  a  whole  frequency  octave  and  should 
thus  be  suitable  for  frequency-comb  generation.  An  in¬ 
teresting  feature  is  that  the  prominent  blue  and  red  spec¬ 
tral  bands  that  are  roughly  mirror  images  of  each  other. 
This  creates  essentially  a  “two-color”  continuum,  and  the 
pulse  in  the  real-space  domain  is  essentially  a  superposi¬ 
tion  of  two  spectrally  wide  sub-pulses  with  their  central 
wavelengths  split  from  the  pump  wavelength  by  a  few 
hundreds  of  nm.  Further  propagation  will  thus  result  in 
spatial  separation  and  a  pair  of  mutually  coherent,  spa¬ 
tially  separated  pulses  will  arise. 


In  conclusion,  we  have  presented  numerical  simulation 
results  simulating  SC  generation  in  the  visible  range  us¬ 
ing  sub-micron  tapered  silica  fibers  which  exhibit  two 
ZDWs.  Our  numerical  simulations  indicate  that  such 
thin  silica  strands,  though  currently  difficult  to  manufac¬ 
ture,  offer  improvements  in  the  quality  of  SC  generated 
spectra,  and  they  should  be  even  more  attractive  white- 
light  sources  than  the  thicker  tapered  fibers  with  a  single 
ZDW  in  the  visible  region. 
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Conclusions 


Long  range  pulse  propagation  in  air 


O) 

c 

0 

13 

o 

o 


0 

0 


0 


O) 

c 

-0 


o 

o 


c 
0 
0  0 

£ 

0  ^3 

-Q  E 

0  <g 
o  J? 

c  Q. 

ro  c 

00  0 


0 

0 


0 

N 

■  MM 

c 

o 

■  HMM 

o 

* 

0 

■O 

0 

E 

0 

0 


o 

0 

o 

E 

c 

0 

D)  0 

°  S. 

£  <D 

0  0 
C  w 
0 


O) 

o 


0 


CL  C  -Q 


H.  Wille  et  al.,  Eur.  Phys.  J.  Appl.  Phys.  20  (2002)  183;  www.teramobile.org 


Control  of  a  chemical  reaction 


Photodissociation  of  CpFe(CO)2CI  (dicarbonylchloro(Af-cyclopentadienyl)iron) 


Ionization  mechanisms 
at  different  field  strengths 


intensity  or  field  strength 


One-electron  atom  in  a  strong  laser  field 


But  no  solution  known  for  the  combination  of  both  fields  !! 


Energy  spectra: 

numerical  simulations  vs.  experiment 


Orbital  symmetry  and  axis  orientation 
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Polyatomic  molecules  ? 
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Ionization  in  intense  laser  field 
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Theoretical  Model 
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•  Field  influence  and  target  properties  are  quasi  separated 


Atomic  ion  signals:  S-matrix  vs.  experiment 
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Angular  distributions  (it-  vs.  a-orbital) 


Node  along  polarization  direction  for 

aligned  molecules  having  HOMO  of  u-symmetry 
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Signature  of  alignment  (7T-HOMO) 
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Dependence  of  total  rates  on  orientation 
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Orientational  dependence  increases  with  intensity 
(consistent  with  highly  nonlinear  nature  of  the  process) 


Dependence  of  total  yields  on  orientation 

No,  800  nm,  45  fs,  2^0*  W/cm2 
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